module eqs

    use polynomial
    use basic
    use spectra
    implicit none

contains

    subroutine a2vd(a,v,d,n,dt,v0,d0)

        real(8), intent(in) :: dt,v0,d0
        integer(4), intent(in) :: n
        real(8), intent(in) :: a(n)
        real(8), intent(out) :: v(n), d(n)

!f2py   intent(in) :: a,n,dt,v0,d0
!f2py   intent(out) :: v,d

        call acc2vd(a,v,d,n,dt,v0,d0)

        return
    end subroutine a2vd

    subroutine targetDispCorr(a,td,ac,n,tp,ntp,ph,pl,dt,v0,d0)

        integer(4), intent(in) :: n
        real(8), intent(in) :: dt,v0,d0
        real(8), intent(in) :: a(n), td(n)
        real(8), intent(out) :: ac(n)
        integer(4), intent(in) :: ntp,ph,pl
        integer(4), intent(in) :: tp(ntp)

!f2py   intent(in) :: n
!f2py   intent(in) :: dt,v0,d0
!f2py   intent(in) :: ntp,ph,pl
!f2py   intent(in) :: a(n), td(n), tp(ntp)
!f2py   intent(out) :: ac(n)

        integer :: i,j,k,pn
        real(8), allocatable :: M(:,:), b(:), tc(:)
        real(8), allocatable :: v(:), d(:), t(:)

        pn = ph - pl + 1
        if ( pn>ntp ) pn = ntp

        allocate(M(ntp,pn))
        allocate(b(ntp))
        allocate(tc(ntp))
        allocate(t(n))
        allocate(v(n))
        allocate(d(n))

        t = [ (((i-1)*dt),i=1,n) ]
        call a2vd(a,v,d,n,dt,v0,d0)

        do k = 1, ntp, 1
            j = tp(k)
            tc(k) = t(j)
            b(k)  = td(j)-d(j)
            do i = pl, pl+pn-1, 1
                M(k,i-pl+1) = dpower(tc(k),i,0)
            end do
        end do
        
        call leastsqs(M,b,ntp,pn)

        ac = a

        do k = 1, n, 1
            do i = pl, pl+pn-1, 1
                ac(k) = ac(k)+b(i-pl+1)*dpower(t(k),i,2)
            end do
        end do

        deallocate(M)
        deallocate(b)
        deallocate(tc)
        deallocate(t)
        deallocate(v)
        deallocate(d)

        return

    end subroutine targetDispCorr

    subroutine detrend(a,ac,n,oh,ol)
        integer(4), intent(in) :: n
        real(8), intent(in) :: a(n)
        real(8), intent(out) :: ac(n)
        integer(4), intent(in) :: oh,ol

!f2py   intent(in) :: n
!f2py   intent(in) :: oh,ol
!f2py   intent(in) :: a(n)
!f2py   intent(out) :: ac(n)

        ac = a
        call polydetrend(ac,n,oh,ol)

        return

    end subroutine detrend

    subroutine basecorr(a,ac,n,oh,ol,dt,v0,d0)
        integer, intent(in) :: n
        real(8), intent(in) :: dt,v0,d0
        real(8), intent(in) :: a(n)
        real(8), intent(out) :: ac(n)
        integer, intent(in) :: oh,ol

!f2py   intent(in) :: n
!f2py   intent(in) :: dt,v0,d0
!f2py   intent(in) :: oh,ol
!f2py   intent(in) :: a(n)
!f2py   intent(out) :: ac(n)
        
        ac = a
        call polyblc(ac,n,oh,ol,dt,v0,d0)

        return
    end subroutine basecorr

    subroutine spectrum(acc,n,dt,zeta,P,np,SPA,INDSPA,SM)
        integer, intent(in) :: n, np, SM
        real(8), intent(in) :: dt,zeta
        real(8), intent(in) :: acc(n), P(np)
        real(8), intent(out) :: SPA(np)
        integer, intent(out) :: INDSPA(np)

!f2py   intent(in) :: n, np, SM
!f2py   intent(in) :: dt,zeta
!f2py   intent(in) :: acc(n), P(np)
!f2py   intent(out) :: SPA(np)
!f2py   intent(out) :: INDSPA(np)

        select case (SM)
            case (1)
                call spafreq(acc,n,dt,zeta,P,nP,SPA,INDSPA)
            case (2)
                call spanmk(acc,n,dt,zeta,P,nP,SPA,INDSPA)
            case (3)
                call spamixed(acc,n,dt,zeta,P,nP,SPA,INDSPA)
            case (4)
                call pspafreq(acc,n,dt,zeta,P,nP,SPA,INDSPA)
            case (5)
                call pspanmk(acc,n,dt,zeta,P,nP,SPA,INDSPA)
            case (6)
                call pspamixed(acc,n,dt,zeta,P,nP,SPA,INDSPA)
            case default
                call spamixed(acc,n,dt,zeta,P,nP,SPA,INDSPA)
        end select

!         SPA = abs(SPA)

        return
        
    end subroutine spectrum

    subroutine spectrumavd(acc,n,dt,zeta,P,np,SPA,INDSPA,SPV,SPD,SPEV,SM)
        integer, intent(in) :: n, np, SM
        real(8), intent(in) :: dt,zeta
        real(8), intent(in) :: acc(n), P(np)
        real(8), intent(out) :: SPA(np), SPV(np), SPD(np), SPEV(np)
        integer, intent(out) :: INDSPA(np)

!f2py   intent(in) :: n, np, SM
!f2py   intent(in) :: dt,zeta
!f2py   intent(in) :: acc(n), P(np)
!f2py   intent(out) :: SPA(np), SPV(np), SPD(np), SPEV(np)
!f2py   intent(out) :: INDSPA(np)

        select case (SM)
            case (1)
                call spavdfreq(acc,n,dt,zeta,P,nP,SPA,INDSPA,SPV,SPD,SPEV)
            case (2)
                call spavdnmk(acc,n,dt,zeta,P,nP,SPA,INDSPA,SPV,SPD,SPEV)
            case (3)
                call spavdmixed(acc,n,dt,zeta,P,nP,SPA,INDSPA,SPV,SPD,SPEV)
            case default
                call spavdmixed(acc,n,dt,zeta,P,nP,SPA,INDSPA,SPV,SPD,SPEV)
        end select

!         SPA = abs(SPA)
!         SPV = abs(SPV)
!         SPD = abs(SPD)
!         SPEV = abs(SPEV)

        return
        
    end subroutine spectrumavd

!     subroutine response(acc,n,dt,zeta,P,ra,rv,rd,SM)
!         integer, intent(in) :: n, SM
!         real(8), intent(in) :: dt, zeta, P
!         real(8), intent(in) :: acc(n)
!         real(8), intent(out) :: ra(n), rv(n), rd(n)

! ! f2py  intent(in) :: n, SM
! ! f2py  intent(in) :: dt, zeta, P
! ! f2py  intent(in) :: acc(n)
! ! f2py  intent(out) :: ra(n), rv(n), rd(n)

!         select case (SM)
!             case (1)
!                 call rfreq(acc,n,dt,zeta,P,ra,rv,rd)
!             case (2)
!                 call rnmk(acc,n,dt,zeta,P,ra,rv,rd)
!             case (3)
!                 call rnmk(acc,n,dt,zeta,P,ra,rv,rd)
!             case default
!                 call rnmk(acc,n,dt,zeta,P,ra,rv,rd)
!         end select
        
!     end subroutine response

    subroutine fdfitspectrum(acc,n,dt,zeta,P,nP,SPAT,a,tol,mit)

        integer, intent(in) :: n, nP, mit
        real(8), intent(in) :: dt, zeta, tol
        real(8), intent(in) :: acc(n), P(nP), SPAT(nP)
        real(8), intent(out) :: a(n)

!p2fy   intent(in) :: n, nP, mit
!p2fy   intent(in) :: dt, zeta, tol
!p2fy   intent(in) :: acc(n), P(nP), SPAT(nP)
!p2fy   intent(out) :: a(n)

        real(8), allocatable :: EP(:), ESPAT(:)

        allocate(EP(nP+2))
        allocate(ESPAT(nP+2))

        EP(2:nP+1) = P
        EP(1) = P(1)*0.5
        EP(nP+2) = P(nP)*1.5
        ESPAT(2:nP+1) = SPAT
        ESPAT(1) = SPAT(1)-(SPAT(2)-SPAT(1))/(P(2)-P(1))*P(1)*0.5d0
        ESPAT(nP+2) = SPAT(nP)+(SPAT(nP)-SPAT(nP-1))/(P(nP)-P(nP-1))*P(nP)*0.5d0

        call fitspectra(acc,n,dt,zeta,EP,nP+2,ESPAT,a,tol,mit)

        deallocate(EP)
        deallocate(ESPAT)
      
    end subroutine fdfitspectrum

    subroutine tdfitspectrum(acc,n,dt,zeta,P,nP,SPAT,a,tol,mit)
        integer, intent(in) :: n, nP, mit
        real(8), intent(in) :: dt, zeta, tol
        real(8), intent(in) :: acc(n), P(nP), SPAT(nP)
        real(8), intent(out) :: a(n)

!p2fy   intent(in) :: n, nP, mit
!p2fy   intent(in) :: dt, zeta, tol
!p2fy   intent(in) :: acc(n), P(nP), SPAT(nP)
!p2fy   intent(out) :: a(n)

        call adjustspectra(acc,n,dt,zeta,P,nP,SPAT,a,tol,mit)
        
    end subroutine tdfitspectrum

end module eqs